clear all;
clc;

% Hochrechnung auf eine Sekunde Bestrahlung (Soll: 10^16 Photonen)
simfaktor = 1e11;
%simfaktor = 1;

% Skala der Colorbar [Gy]
gray = 1e-2;

if(exist('filename', 'var') == 0)
    [filename, pathname] = uigetfile('.log');
    file = fopen(fullfile(pathname, filename), 'r');
    
    % Name xCoord yCoord zCoord eDeposit
    input = fscanf(file,'%*s\t%e\t%e\t%e\t%e\n',[4 inf]);
    
    % Bitte daran denken das File zu schliessen, sonst gibt es Probleme
    fclose(file);
    % Eventuell an diese Stelle ein 'close all'
end

% Anzahl Zellen
cells = 50;

% Dichten der Teilnehmer [kg/m^3]
rho_tumor       = 1200;
rho_brain       = 1000;
rho_eyes        = 1000;
rho_afterbrain  = 1000;
rho_muscle      = 1000;
rho_skin        = 1000;
rho_bone        = 1000;

% Temporär: Energie in Joule
input(4,:) = input(4,:).*1.602e-16;

% Input auf Zellengroesse diskretisieren
input(1:3,:) = floor(input(1:3,:) .* (cells-1));

% Maximum und Minimum der x-, y- und  z-Werte finden
minimum_x = min(input(1,:));
maximum_x = max(input(1,:)) - minimum_x;

minimum_y = min(input(2,:));
maximum_y = max(input(2,:)) - minimum_y;

minimum_z = min(input(3,:));
maximum_z = max(input(3,:)) - minimum_z;

% Input normieren
input(1,:) = input(1,:) - minimum_x;
input(2,:) = input(2,:) - minimum_y;
input(3,:) = input(3,:) - minimum_z;


% Schrittgroessen finden
dx = maximum_x/(cells - 1);
dy = maximum_y/(cells - 1);
dz = maximum_z/(cells - 1);

% Deposition in Gray
input(4,:) = input(4,:)./(rho_tumor.*dx.*dy.*dz.*0.001^3);

% Film festlegen mit 10fps
film = avifile('radiation_dosimetry.avi', 'quality', 100, 'fps', 10);

% Figure Properties setzen
fig = figure('Position', [0, 0, 800, 800]);
load cmap;
colormap(cmap);

% Fuer jede Zelle
for i=1:cells
    matrix = zeros(cells, cells);
    % Fuer jede Inputzeile
    for j=1:size(input,2)
        if(input(2,j) >= (i-1)*dy && input(2,j) <= i*dy)
            matrix(ceil(input(1, j)/dx + 0.5), ceil(input(3, j)/dz + 0.5)) = matrix(ceil(input(1, j)/dx + 0.5),...
                ceil(input(3, j)/dz + 0.5)) + input(4, j);
        end
    end
    
    % Falschfarbenplot
    surface_handle = pcolor(gca, matrix);
    set(surface_handle, 'EdgeColor','none');
    title('Radiation Dosimetry of Patient');
    ylabel('zAxis')
    xlabel('xAxis')
    set(gca, 'CLim', [0, gray./simfaktor]);
        
    hcb = colorbar;

    axis equal tight;

    %frame = getframe(fig);
    %film = addframe(film, frame);
    print('-dpng', sprintf('output/deposition_%02d',i));
end

film = close(film);
close all;

